Monte Carlo simulation of nonlinear Couette flow in a dilute gas 



Jose Maria MontanercO 
Departamento de Electronic/! e Ingenieria Electromecdnica, 
Universidad de Extremadura, E-06071 Badajoz, Spain 

Andres Santos0 and Vicente GarzoEl 
Departamento de Fisica, Universidad de Extremadura, 
E-06071 Badajoz, Spain 
(February 1, 2008) 

The Direct Simulation Monte Carlo method is applied to solve the Boltzmann equation in the 
steady planar Couette flow for Maxwell molecules and hard spheres. Nonequilibrium boundary con- 
ditions based on the solution of the Bhatnagar-Gross-Krook (BGK) model for the Couette flow are 
employed to diminish the influence of finite-size effects. Non-Newtonian properties are characterized 
by five independent generalized transport coefficients: a viscosity function, a thermal conductivity 
function, two viscometric functions, and a cross coefficient measuring the heat flux orthogonal to the 
thermal gradient. These coefficients depend nonlinearly on the shear rate. The simulation results 
are compared with theoretical predictions given by the Grad method and the BGK and the ellip- 
soidal statistical (ES) models. It is found that the kinetic models present a good agreement with 
the simulation, especially in the case of the ES model, while the Grad method is only qualitatively 
reliable for the momentum transport. In addition, the velocity distribution function is also measured 
and compared with the BGK and ES distributions. 

PACS numbers: 47.50.+d, 51.10.+y, 05.20.Dd, 05.60.-k 



I. INTRODUCTION 

One of the most interesting states for analyzing transport processes far from equilibrium is the steady planar Couette 
flow. The physical situation corresponds to a system enclosed between two infinite, parallel plates in relative motion 
and, in general, kept at different temperatures. These boundary conditions lead to combined heat and momentum 
transport. If x and y denote the coordinates parallel to the flow and orthogonal to the plates, respectively, then the 
corresponding steady hydrodynamic balance equations are 
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where u — u x x is the flow velocity, P is the pressure tensor, and q = q x x + q y y is the heat flux. The presence of q y in 
Eq. (Q) indicates that a thermal gradient BT/8y is induced by the velocity gradient, even if both plates are kept at 
the same temperature. The balance equations (0) and (||) do not constitute a closed set unless the dependence of the 
pressure tensor and the heat flux on the hydrodynamic fields is known. If the gradients are small, the fluxes P and q 
are described by the Navier-Stokes (NS) constitutive relations, which in this problem yield 
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q x = 0, q y = ~ K 0-g-, (4) 

where t/q and kq are the NS shear viscosity and thermal conductivity coefficients, respectively. As a consequence of 
the absence of normal stress differences in the NS description, the hydrostatic pressure p = (P xx + P yy + Pzz)/3 is a 
constant, on account of the balance equation (|j). 

Even in the linear regime described by the NS equations, one still needs to know the spatial dependence of the 
transport coefficients to obtain the exact solution of the hydrodynamic equations. The problem becomes tractable 
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in the case of a low density gas, where the state of the system is completely specified by the velocity distribution 
function f(r,v;t), which obeys the Boltzmann equations A relevant dimensionless quantity in a dilute gas is the 
Knudsen number Kn = \/ih, defined as the ratio of the mean free path A to the scale length of the hydrodynamic 
gradients 1%. In many laboratory conditions, Kn -C 1 and so the Boltzmann equation can be solved by means of the 
Chapman-Enskog method as an expansion of the distribution function in powers of the Knudsen number The zeroth 
order approximation leads to the Euler hydrodynamic equations, while the first order approximation yields the NS 
equations with explicit expressions for the transport coefficients rjo and kq. The results show that the ratio tjq/ko is 
a constant. Consequently, it then follows from the NS hydrodynamic equations (|l])-(f§) that the flow velocity profile 
is quasi-linear, 

Vo-^ = const, (5) 
dy 

and the temperature is quasi-parabolic, 

Kq—\ T = -—[rio—^] = const. (6) 

Note that the profile of u x is not strictly linear, due to the space dependence of r/o through the temperature. Anal- 
ogously, the temperature profile is not strictly quadratic. In fact, the specific form of both profiles depends on the 
interaction potential under consideration. On the other hand, from Eqs. (||) and (^) it is easy to derive a nice result, 
namely that if the temperature T is seen as a function of u x rather than as a function of the coordinate space y, then 
one has 

— = (7) 

dul K 

This is a sort of nonequilibrium "equation of state," according to which the temperature is a quadratic function of 
the flow velocity. Moreover, the "curvature" of the profile is practically universal, given the weak influence of the 
interaction potential on the Prandtl number Pr = 5fcs?7o/2mKo — §> where fee is the Boltzmann constant and m is 
the mass of a particle. 

Since the mean free path is inversely proportional to the density, the Knudsen number, at a given value of the scale 
length £h, increases as the gas becomes more rarefied. In general, when the Knudsen number is not small, the NS 
relations are not expected to hold and the transport must be described by nonlinear constitutive equations. In the 
special case of Maxwell molecules (particles interacting via an r~ 4 potential) , it has been showroo that the Boltzmann 
equation admits a consistent solution in the nonlinear Couette flow characterized by a constant pressure p and profiles 
similar to those obtained in the NS regime, Eqs. (Eh - ©) except that 770 and kq are replaced by a generalized shear 
viscosity coefficient 77(a) = rjaF ri (a) and a generalized thermal conductivity coefficient n yy (a) — kqF k {o), respectively. 
Here, a = (r)o/p)du x /dy is a constant (dimensionless) shear rate and F v and F K are nonlinear functions of a. In 
addition, P xx ^ P yy P zz and q x 7^ 0. In this problem, the hydrodynamic scale length can be identified as 
£h ~ \J ksT I m(du x I dy)~ x , while the mean free path is A ~ ksT / m{rjQ / p) . Thus, the reduced shear rate a 
represents the Knudsen number in this problem, i.e. a ~ Kn. Henceforth, we will use the reduced shear rate a to 
refer to the Knudsen number Kn. The solution considered in Refs. ||,^] describes heat and momentum transport for 
arbitrary velocity and thermal gradients in the bulk domain, where boundary effects are negligible. On the other 
hand, the full nonlinear dependence of F n (a) and F K (a) is not explicitly known, since it involves the infinite hierarchy 
of moment equations. Their knowledge is limited to super-Burnett order and the result isu F r) (a) = 1 — 3.111a 2 and 
F K (a) = 1 - 7.259a 2 . 

Consequently, if one wants to get the transport properties for arbitrary values of the shear rate and the thermal 
gradient, one must resort to approximate schemes or to computer simulations. In the first alternative, explicit 
expressions for the nonlinear transport coefficients in the Couette, flow have been obtained from exact solutions of 
the Bhatajagar-Gross-Krook (BGK) modelot] and related modelsLTEl for general interactions, as well as from the Grad 
method.E2l In the simulation side, Risso and Corderdl3 have recently studied the shear-rate dependence of F v and F K 
by means of molecular dynamics simulations of a hard disk gas. Comparison between the different analytical results 
with those obtained from the simulation shows that the predictions given by kinetic models, are in better agreement 
than those given by the Grad method, especially in the case of the thermal conduct ivity.u Nevertheless, given the 
difficulties associated with molecular dynamics simulations to achieve large shear rates in a dilute gas, the above 
comparison is restricted to a range of shear rates for which non-Newtonian effects are hardly significant. For instance, 
the shear viscosity has only decreased around-40% with respect to its Navier-Stokes value for the largest value of the 
shear rate considered by Risso and Cordero.t2l In order to overcome such difficulties and extend the range of values 
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of a, one may use the so-called Direct Simulation Monte Carlo (DSMC) method,0 which is known to qualify as an 
efficient and accurate method to numerically solve the Boltzmann equation. 

The aim of this paper is to solve the Boltzmann equation by means of the DSMC method for a gas subjected to 
the planar Couette flow. The motivation of this work is twofold. On the one hand, we want to test the reliability of 
the far from equilibrium results obtained from kinetic models and the Grad method by making a comparison with 
the Boltzmann solution in the case of Maxwell molecules, for which the form of the hydrodynamic profiles in the 
bulk region (far from the boundaries) is known. We will determine not only the hydrodynamic profiles but also the 
nonlinear transport coefficients and the velocity distribution function. On the other hand, we want to investigate 
whether the above results for a system of Maxwell molecules extend to other interaction potentials. This extension 
holds when the Boltzmann equation is replaced by kinetic models where, in terms of a conveniently scaled space 
variable, all the results are independent of the interaction law. Thus, we will also solve numerically the Boltzmann 
equation by the DSMC method for a hard-sphere gas. 

Since we are interested in describing transport properties in the bulk region, i.e., free from finite-size effects, we need 
to use appropriate boundary conditions in the simulations. In the conventional boundary conditions ,113113 the gas is 
assumed to be enclosed between two baths at equilibrium in relative motion and, in general, at different temperatures. 
Under these conditions, a particle leaving the system is formally replaced by a particle coming from the bath, so 
the incoming velocity is sampled from a local equilibrium distribution. As a consequence, there exists a mismatch 
between the velocity distribution function of the reemitted particles and that of those particles of the gas located near 
the wall and moving along the same direction. In this case, in order to inhibit the influence of boundary effects one 
needs to take very large systems (normal distance between the plates much larger than the mean free path), what 
is not practical from a computational point of view. To overcome this difficulty, a possibility is to assume that both 
baths are out of equilibrium in a state close to that of the actual gas. Since such a state is not known "a priori"™ 
in this paper we assume that the state of the baths is described by the BGK solution of the planar Couette flow.H 
Although the boundary effects are still unavoidable, we expect that the above mismatch between reemitted and gas 
particles will be much smaller. As a matter of fact, the use of these conditions has been shown-tp be appropriate to 
analyze bulk transport properties in the special case of planar Fourier flow (both walls at rest). 113 

The plan of the paper is as follows. In Sec. [n| we give a brief description of the planar Couette flow and a summary 
of the main results obtained from the Boltzmann equation and kinetic models. The boundary conditions used in the 



simulations and the DSMC method are described in Sec. IB. Section IV presents the main results of the paper, where 



special attention is paid to the nonlinear transport coefficients. A comparison with the analytical results derived from 
the BGK and the ellipsoidal statistical (ES) models and from the Grad method is also carried out. The comparison 
shows in general a good agreement of the kinetic models with computer simulations, even for large shear rates. In 
addition, the velocity distribution function obtained from the simulation in the bulk domain is compared with the 
ones given by the kinetic models. It is shown again that the agreement is qualitatively good. We close the paper in 
Sec. with some concluding remarks. 



II. DESCRIPTION OF THE PROBLEM AND SUMMARY OF THEORETICAL RESULTS 

Let us consider a dilute gas. In this case, a kinetic description is sufficient to characterize the state of the system 
by means of the velocity distribution function /(r, v;t). This distribution function obeys the nonlinear Boltzmann 
equation, which in the absence of external forces is given byu 

^+v-V/=J[/,/], (8) 

where 

J[fJ] =J d viJ dkgl(g,k) [f(v')f(v[) - /(v)/( Vl )] (9) 

is the collision operator. In this equation, I(g. k) is the differential cross section, g = |v— vi | being the relative velocity, 
and (v'jV^) are precollisional velocities yielding postcollisional velocities (v, Vi). From the distribution function, one 
may define the hydrodynamic quantities 

n= fdvf, (10) 

u=-fdvvf, (II) 
n I 
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~nfc B T=| jaVV 2 /, (12) 

and the momentum and heat fluxes 

'dvW/, (13) 



q=- J dvV'Vf. (14) 

Here, n is the number density, u is the flow velocity, T is the temperature, P is the pressure tensor, q is the heat flux, 
and V = v u is the peculiar velocity. In addition, the equation of state is that-pf an ideal gas, i.e., p = nksT. 

Most of the known solutions to Eq. (||) for spatially inhomogeneous statesli-3 correspond to the special case of 
Maxwell molecules, namely, a repulsive potential of the form V(r) ~ r~ 4 . For this potential, the collision rate gl(g, k) 
is independent of the relative velocity and this allows the infinite hierarchy of velocity moments to be recursively 
solved in some specific situations. Furthermore, the NS transport coefficients 770 and kq can be exactly obtained from 
the Chapman-Enskog method.cl They are given by 

P 15fci3 

Vo = ~, k = rfo, (15) 

v 4 m 

where v — On, 6 being an eigenvalue of the linearized Boltzmann collision operator0 

In this paper we are interested in studying the planar Couette flow for a dilute gas. We consider a gas enclosed 
between two parallel plates in relative motion and maintained at different temperatures. Under these conditions, the 
system is driven out of equilibrium by the combined action of the velocity and thermal gradients along the direction 
normal to the plates. After a transient period, the gas is expected to reach a steady state and the corresponding 
Boltzmann equation reads 

Vy-^f = J[fJ], (16) 

where we have chosen the axis y as the one normal to the plates. In general, this equation must be salved subjected 
to specific boundary conditions. Nevertheless, in the same spirit as in the Chapman-Enskog method,EI one may look 
for "normal" solutions in which all the space dependence of the distribution function occurs through a functional 
dependence on the hydrodynamic fields. The normal solution describes the state of the gas in the hydrodynamic 
regime, namely, for times much longer than the mean free time and for distances from the walls much larger than 
the mean free path. As mentioned in the Introduction, it has been provedaa that, in the case of Maxwell molecules, 
Eq. (|l6|) admits an exact solution corresponding to the planar Couette flow problem. This solution belongs to the 
normal class, so that no explicit boundary conditions appear. In other words, all the space dependence of / is given 
through the local density, the local velocity, the local temperature, and their gradients. The solution is characterized 
by hydrodynamic profiles that are a simple generalization of those predicted by the NS approximation, Eqs. (^) _ (0), 
namely 

p = const, (17) 
1 dur. 



v(y) dy 
1 d 



v(y) dy 



a = const, (18) 
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T=-Pr— 7 (a). (19) 



The dimensionless parameter 7(a) is a nonlinear function of the reduced shear rate a that, by construction, behaves 
as 7 w a 2 /5 in the limit a — > 0. Again, the temperature can be seen as a quadratic function of the flow velocity, 
but now the coefficient ?7o/ K o appearing in Eq. (Q) is replaced by a shear-rate dependent coefficient (rj Q /K Q )5j(a)/a 2 . 
Furthermore, in this solution the pressure tensor is independent of the thermal gradient and the heat flux vector is 
linear in the thermal gradient, but all these fluxes are nonlinear functions of the shear rate. This nonlinear dependence 
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can be characterized through five generalized transport coefficients. First, the shear stress P xy defines a generalized 
shear viscosity 77(a) as 



du x du x 
— ^- V0 F v (a) — 



P xv = -ri{a)-^ = -r l0 Fr l {a)-^. (20) 



Analogously, the component of the heat flux parallel to the thermal gradient defines a generalized thermal conductivity 
coefficient n yy (a): 

dT dT 
% = - K yy( a )-Qj; = -KoF^a) — . (21) 

The dimensionless functions F ri (a) and F K (a) are the most relevant quantities of the problem. They are related to 
the function 7(a) by 7(a) = a 2 F JI (a)/5F K (a). Normal stress differences are different from zero and are measured by 
two viscometric functions ^1.2(0) defined by 

'' J:l = ^(a)a 2 , (22) 



Pzz Pyv = * 2 (a)a 2 . (23) 
V 

Another interesting quantity related to the pressure tensor is the friction coefficient /_i(a) = —P X y/P yy . To NS order, 
we simply have /i(a) = a. In the non-Newtonian regime, we generalize this coefficient as /z(a) = aF^(a), where the 
friction function F^^a) is 

F -<"> = i-pJ-'Uiw < 24 > 

Finally, there exists a nonzero component of the heat flux orthogonal to the thermal gradient given by 

d d 

1x = -Kx V (a)—T = -K $(a)a— T. (25) 
oy ay 

The three functions ^ 1,2(0) and $(a) are generalizations of Burnett coefficients. In fact, \I/i(0) ^ —14/5, ^(O) = 4/5, 
and ^(O) = —7/2 for Maxwell molecules and for hard spheres in the first Sonine approximation.^ The determination of 
the nonlinear transport coefficients F n , F K , ^i^, and would imply the solution of an infinite hierarchy that cannot 
be solved in a recursive way.Q This hierarchy can only be salved step by step when one performs a perturbation 
expansion in powers of the shear rate. In fact, Tij and Santoso determined the solution up to super-Burnett order: 

F n (a) = 1 - 3.111a 2 + 0(a 4 ), (26) 



F K (a) = 1 - 7.259a 2 + C(a 4 ). (27) 



FromEqs. (|26|) and (|27|) it follows that 7(a) = (a 2 /5)[l + 4.148a 2 + C(a 4 )]. In addition, F^a) = 1 - 1.911a 2 + C(a 3 ). 

Although the above analyses are valuable, they have two main limitations. On the one hand, they are restricted to 
the special case of Maxwell molecules. For other interaction potentials (e.g., hard spheres), the collisional moments 
involve all the moments of the distribution function and, as a consequence, the hydrodynamic profiles (p^7|)- (|l9|) are 
not strictly true. On the other hand, even for Maxwell molecules, the above perturbative solution is not useful 
for finite shear rates. These two limitations can be overcome, in the context of analytical methods, by introdij 
additional approximations, such as the Grad method,t3 or by describing the system by means of kinetic models! 
In these approaches, one looks for a solution having the same hydrodynamic profiles as in the case of the Boltzmann 
equation, cf. Eqs. (p^7|)-(^9|). As before, this solution describes the properties of the system in the bulk region, which 
is insensitive to the details of the boundary conditions. From the Grad method and from the BGK and ES kinetic 
models, one explicitly gets the full nonlinear shear-rate dependence of the transport coefficients in a non-perturbative 
way. Moreover, the results are universal, namely the functions F v (a), F K (a), . . .are independent of the interaction 
potential, provided that the reduced shear rate is defined as in Eq. (|18[) with v = p/rjo- 
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(28) 



(29) 



is the local equilibrium distribution function. When the approximation (|28|) is inserted into the Boltzmann equation 
(||) and velocity moments are taken, one gets a closed set of equations for n, u, P, and q. According to the geometry 
of the planar Couette flow, the re ., ar e eight independent moments instead of the original thirteen moments appearing 
in Eq. (^8|). Risso and CordercEJiij found that the set of independent moment equations, neglecting nonlinear terms 
in the fluxes, admits a solution consistent with the profiles (pL T|) — ( 19). In addition, they obtained explicit expressions 
for the transport coefficients introduced above as nonlinear functions of the reduced shear rate a. Those expressions 
are displayed in the Appendix. 

Now we consider the kinetic model approach. The basic idea is to replace the detailed Boltzmann collision operator 
by a much simpler term that otherwise retains the main physical properties of the original operator J[f, /]. The most 
familiar choice is the BGK model,EJ 



4f,/]--K/-M 



(30) 



where v is an effective collision frequency that depends on the temperature, according to the interaction potential. 
The NS transport coefficients of the BGK model are ijq = p/v and k,q = 5pk B /2mv. Thus the BGK model has the 
drawback that it predicts an incorrect value for the Prandtl number, Pr = 1. This is a consequence of the fact that v 
is the wly adjustable parameter in the model. This deficiency is corrected by the so-called ellipsoidal statistical (ES) 
model,B in which case 
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where £ is again an effective collision frequency and the reference function fa is 

fa = mr- 3/2 (det A)~ 1/2 exp (-A" 1 : VV) , 



(31) 



(32) 



where Aij = (2k B T/m)Pr~ 1 5ij — 2(Pr _1 — l)Pij/mn. The NS coefficients are now ijq =p/((Pr~ 1 ) and kq = 5pk B /2m£. 
If, as in Eq. ([To]), we define a collision frequency v = p/rjo, then v = (Pr -1 in the ES model. Note that the ES model 
reduces to the BGK model if we formally make Pr = 1. Therefore, the ES model can be seen as an extension of the 
BGK model to account for the correct Prandtl number, which can be particularly relevant in the Couette flow where 
heat transport and momentum transport are coupled. The results derived from the BGK and the ES models for the 
Couette flow problem are also given in the Appendix. Apart from obtaining the transport properties, the velocity 
distribution function can be explicitly written. In particular, the BGK distribution is given by 
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Here, (i ,*i) = (0, 1) if £ y > and (t , h) = [1,2/(1 - a)} if £ y < 0. Besides, £ = (m/2k B T) 1 / 2 V , 
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(34) 



(35) 



is a (local) reduced thermal gradient. Equation ( |33j ) shows that the distribution function is a highly nonlinear function 
of the reduced gradients a and e. 
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In Ref . |g a comparison between the analytical results obtained from the Grad method and the BGK and ES models 
with those obtained from molecular dynamics simulationsEJ for hard disks was carried out. It was found that the 
three theories reproduced quite well the simulation data for F v , but the Grad method failed for F K and $. The 
latter quantity was reproduced better by the ES model than by the BGK model. Notwithstanding this, more definite 
conclusions could not be reached because the simulation data were restricted to rather small shear rates, namely 
a < 0.2. In this range of shear rates, non-Newtonian effects are not especially significant. In addition, although the 
molecular dynamics simulations correspond to a very dilute gas (area fraction <f> « 1%), the collisional contributions 
to the transport coefficients (absent in a Boltzmann description) are not strictly zero. Finally, as a minor point, 
conclusions drawn in the context of two-dimensional systems should not be extrapolated without caution to the more 
realistic case of three dimensions. As said in the Introduction, the aim of this paper is to solve numerically the 
Boltzmann equation by means of the DSMC method for the planar Couette flow and compare the results with the 
theoretical predictions. We will consider three-dimensional systems of Maxwell molecules and hard spheres subjected 
to shear rates as large as a ~ 1.2. In addition to the nonlinear transport coefficients, the comparison will be extended 
to the level of the velocity distribution function itself. 



III. BOUNDARY CONDITIONS AND MONTE CARLO SIMULATION 



A. Boundary conditions 

The goal now is to solve numerically the Boltzmann equation corresponding to the planar Couette flow by using 
the successful DSMC method.El The gas is enclosed between two parallel plates located at y = and y = L, which 
are moving along the x-direction with velocities Uo = Uqx and Ui = C/lx, respectively. In addition, they are kept 
at temperatures To and Tj,, respectively. In order to solve Eq. (|l6|), we need to impose the corresponding boundary 
conditions. They can be expressed in terms of the kernels Ko^(v, v ') defined as follows. When a particle with velocity 
v' hits the wall at y = L, the probability of being reemitted with a velocity v within the-range dv is Kl(v, v')dv; 
the kernel Kq(v, v') represents the same but at y = 0. The boundary conditions are thenO 



e(±v y )\v v \f(y = {0,L},v) = e(±v v ) J dv' \v' y \K , L (v,v') 

xG(Tv' y )f(y={0,L}y,t). (36) 

In this paper we consider boundary conditions of complete accommodation with the walls, so that Kq^(v,v') = 
-Ko,l(v) does not depend on the incoming velocity v' and can be written as 

K 0<L (v) = A-le{±v y )\v y \ct> 0tL (v), A . L = f dvQ(±v y )\v y \^ L (v), (37) 

where <fto.z,(y) represents the probability distribution of a fictitious gas in contact with the system at y = {0, L}. 
Equation ( |37| ) can then be interpreted as meaning that when a particle hits a wall, it is absorbed and then replaced 
by a particle leaving the fictitious bath. Of course, any choice of 4>o,l(v) must be consistent with the imposed wall 
velocities and temperatures, i.e., 

U ,l = dvv x (j) . L (v), (38) 



k B T ,L = -m J dv (v - U o ,l) 2 0o,l(v). (39) 
The simplest and most common choice is that of a Maxwell-Boltzmann (MB) distribution: 

Under these conditions, the system is understood to be enclosed between two independent baths at equilibrium. While 
the conditions (^0|) are adequate for analyzing boundary effects,ErE3 they are not very efficient when one is interested 
in obtaining the transport properties in the bulk region. In order to inhibit the influence of boundary effects, it is 
more convenient to imagine that the two fictitious baths are in nonequilibrium states resembling the state of the actual 



m(v - U ,l) 



2fcsT 0i L 



(40) 
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gas near the walls. More specifically, we can assume that the fictitio us g ases are described by the BGK equation, 
whose exact solution for the steady planar Couette flow is given by Eq. (|33| ) . In this case, the probability distributions 
<?Vl(v) are 



/BGK/ \ —3/2 m «0,Z(1 + "0,i) 3/2 f' 1 

9o,l l v J — n 



£o,L\ v y\ 



^ dt [2t- (1 - a QtL )t 2 \ 

•I to 



-5/2 



x exp ■ 



2k B T r 



1/2 



2a ,. 



1 - t 



1 + a 0t L 



2a a , L 1 \ ■:> > 

Wrr - Uq, L + — ) + V„ + V 

1 + a , L £o,. 



1 + olq.l £o,LVy 2k B T 0iL 2t - (1 - a ,L)i 2 

2 



(41) 



where (t , ti) = (0, 1) if v y > and (t , ti) = [1, 2/ (1 - cv ,l)] if v y < 0, and a ,L = £o,l/[£o,l + 87bgk(o')] 1/2 - Here, 
a' is the estimated value of the reduced shear rate, as predicte d b y the BGK model for specific values of the boundary 
parameters Uq,l and To,l, and 7bgk(i') is obtained from Eq. ( A6). Given the values of the four independent boundary 
parameters Uq.l and To.l (as well as the distance L), the shear rate a' and the local thermal gradients £q,l are fixed 
by the conditions (17)-(|l9|). Therefore, 

UL ~ U0 (42) 



1 [2k B T Q , L 
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In these equations, C is related to the actual separation L between the plates through the nonlinear equation 



L = 



ds 



(43) 



(44) 



where s is a variable in terms of which the temperature is a quadratic function, namely T(s) = Tq[1 + 
eo(m/2kBTo) 1 / 2 s — m'yBGK(a')Prs 2 /ksTQ], and the s-dependence of the collision frequency v(s) appears only through 
the temperature (taking into account that p — const). The solution of the nonlinear set of equations (^)-(Q) gives 
a ' , eo.L, and C for any choice of t/o,i, Tb^, and L. Nevertheless, from a practical point of view, it is more convenient 
to fix Uq, T ^ l , a', and e as independent parameters. Without loss of generality we take U = and T L = 1. In 
addition, we will choose eo = 0. This implies that, if boundary effects were absent, the temperature would have a 
maximum at the lower plate y = 0. Thus, 



U L = a'C, C 



A 



2Pr 7B GK(a') 



1/2 



—4 



(45) 



(46) 



In the above equations A = To — 1 and we have taken m = 1 and ks = \- Finally, the actual distance L is given by 
Eq. (44). However, from the simulation point of view, it is more useful to express L in terms of the collision frequency 
v corresponding to the temperature Tl and the average density n. In other words, instead of Eq. (^) we use 

Ti = — [ dyn(y) — — [ ds ' 



V (8) 

to determine L. For the sake of concreteness, let us consider repulsive potentials, for which v 
where u> ranges from (Maxwell molecules) to \ (hard spheres). In that case, 



L = 



1 



ds[T{s)\ 



(47) 

V(n/n)(T/T L r, 
(48) 



For Maxwell molecules, this simply reduces to L — L/v, while for hard spheres one has L — (C/v) tan _1 (vA)/vA- 
In summary, given the values of a' and A, the separation L between the plates and the velocity of the upper plate 
Ul are uniquely determined. In addition, the thermal gradient at the upper plate cl is also determined, while the 
value at the lower plate is fixed as eo = 0. The knowledge of these boundary parameters allows one to obtain the 
distributions ^ K (v) , according to Eq. ( [4l| ) . 
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B. The DSMC method 



Now, wc briefly describe the numerical algorithm we have employed to solve the Boltzmann equation by means of 
the so-called Direct Simulation Monte Carlo (DSMC) method.liil In this method, the velocity distribution function is 
represented by the velocities {v^} and positions {?/,} of a sufficiently large number of particles N . Given the geometry 
of the problem, the physical system is split into layers of width Sy, sufficiently smaller than the mean free path. The 
velocities and coordinates are updated from time t to time t + St, where the time step St is much smaller than the mean 
free time, by applying a streaming step followed by a collision step. In the streaming step, the particles are moved 
ballistically, yi — > yi + v iy St. In addition, those particles crossing the boundaries are reentered with velocities sampled 
from the corresponding probability distribution _Ko,l( v )- Suppose that a particle i crosses the lower plate between 
times t and t + St, i.e., yi + Vi y 5t < 0. Then, regardless of the incoming velocity Vj, a new velocity Vj is assigned 
according to the following rules. First, a velocity Vj (with Vi y > 0) is sampled with a probability proportional to 
|u a |<^ B (v). If one is considering the MB boundary conditions, this velocity is accepted directly. Otherwise, the above 
acts as a "filter" to optimize the acceptance-rejection procedure and the velocity Vj is accepted with a probability 
proportional to the ratio BGK (v)/</>^ B (v). If the velocity is rejected, a new velocity Vj is sampled and the process 
is repeated until acceptance. The new position is assigned as Vi y (St + yi/vi y ). The process is analogous in the case of 
the upper plate. 

The collision step proceeds as follows For each layer a, a pair of potential collision partners i and j are chosen 
at random with equiprobability. The collision between those particles is then accepted with a probability equal to 
the corresponding collision rate times St. For hard spheres, the collision rate is proportional to the relative velocity 
|vj — Vj|, while it is independent of the relative velocity for Maxwell molecules (an angle cut-off is needed in the 
latter case). If the collision is accepted, the scattering direction is randomly chosen according to the interaction law 
and post-collisional velocities are assigned to both particles, according to the conservation of momentum and energy. 
After the collision is processed or if the pair is rejected, the routine moves again to the choice of a new pair until the 
required number of candidate pairs has been taken. 

In the course of the simulations, the following "coarse-grained" local quantities are computed. The number density 
in layer a is 



N a nL X 



n 



(N/L)Sy NSy ^ 



where Q a (y) is the characteristic function of layer a, i.e., Q a (y) = 1 if V belongs to layer a and is zero otherwise. 
Similarly, the flow velocity, the temperature, the pressure tensor, and the heat flux of layer a are 

I N 

Ua = aO/i) v i. ( 50 ) 

a 2 — 1 

N 

k B T a = — = — - Ve Q (y;)(v 2 -u Q ) 2 , (51) 
n a 3N a f-f 

2—1 

L N 

P a = j^my]Q a (yi)(vi - u n )(v, - u a ), (52) 
L m N 

Qa = ^^yE©a(2/i)( v i- U a) 2 ( v i- U a)- (53) 

From these quantities one can get local values of the gradients and of the transport coefficients. For instance, the 
reduced shear rate is 

a a = — [ 7 ±) ^ 54 

vn a \T a J by 

and the viscosity function is 







F Via = -^L. (55) 

As said before, in the simulations we take units such that m = 1, Tjj = 1, and ks — \- It remains to fix the time 
unit or, equivalently, the length unit. The standard definition of mean free path in the case of hard spheres iscl 

\/2nir<7 2 ' ^ 

where a is the diameter of the spheres. The Navier-Stokes shear viscosity is (in the first Sonine approximation) 
rja = 5(TOfcsT/7r) 1 / 2 /16cr 2 . Consequently, the effective collision frequency v — p/r)o and the mean free path A are 
related as v = (8/50r)(2fc B r/™) 1 / 2 /A. As usual, we choose the mean free path corresponding to the average density 
n as the length unit. This in turn implies that V — 8/5^/7? ~ 0.903. For convenience, we take the latter value for 
Maxwell molecules as well. The typical values of the simulation parameters are N — 2 x 10 5 particles, a layer width 
Sy = 0.02, and a time step St = 0.003. 

The procedure to measure the relevant quantities of the problem is as follows. First the values of the imposed shear 
rate a' and the temperature difference A are chosen. This choice fixes the system size L, as well as the upper velocity 
Ul and the upper thermal gradient €l, according to Eqs. (ft5|), jiq), and (|4g|). Starting from an equilibrium initial 
state with T = Tq, the system evolves driven by the boundary conditions described before. After a transient period 
(typically up to t = 25), the system reaches a steady state in which the quantities fluctuate around constant values. 
In this state, the balance equations predict that the quantities u y , P xy , P yy , and u x P xy + q y are spatially uniform and 
this is used in the simulations as a test to make sure that the steady state has been achieved. Once the steady state is 



reached, the local quantities (49)-l55l) are averaged over typically 100 snapshots equally spaced between t = 25 and 



t = 55, which corresponds to about 60 collisions per particle in the case of hard spheres. 

C. Test of the numerical algorithm 

Before closing this Section, it is worthwhile carrying out a test of the reliability of the numerical method. To that 
end, we have simulated the BGK equation by a DSMC-like method similar to the one described in Ref. |l9|. If the 
boundary conditions are implemented correctly and the simulation parameters are well chosen, then the simulation 
results should agree with the theoretical BGK predictions. We have checked that this indeed the case. As an example, 
consider the hard-sphere situation with a' = 1 and A = 5. This corresponds to 7bgk = 0.248, L — 1.81, Ul — 3.17, 
and 6l = —3.15, where in this case Pr = 1. Figure |l| shows the marginal velocity distributions of particles reemitted 
by the walls, 

KoAZv) = ( 2kBTo > L ) 1 r dv x r dv z K , L (-v), (57) 



as functions of £ y = (m/2fcsTo ! L) 1 / 2 u a . The case £ y > (£ y < 0) corresponds to particles that are reemitted from the 
lower (upper) plate. The agreement with the imposed distribution is excellent. Note the strong asymmetry between 
the distributions corresponding to both plates, in contrast to the symmetry of the MB distributions obtained from 
(pp|). The temperature and velocity profiles are shown in Fig. || The simulation values overlap, within statistical 
fluctuations, with the theoretical predictions. Apart from the profiles, we have verified that the generalized transport 
coefficients obtained from the fluxes also agree with the theory. A more stringent test is provided in Fig. ||, where the 
marginal distribution function 

1 /2k B T^ 1/2 z -00 r°° ' - ^ 1/2 



^y) = -[~it) J_Jv x J_Jv z m, €v=^J «, (58) 

is plotted at the point y/L = 0.5, which corresponds to a local thermal gradient e = —0.60. It is apparent again that 
an excellent agreement exists between simulation and theory. 

IV. RESULTS 

This Section is devoted to a comparison between the simulation results for the Boltzmann equation obtained by the 
simulation method described in the previous Section and the theoretical predictions provided by the Grad method 
and the BGK and ES kinetic models. The comparison will be carried out at the levels of the transport coefficients 
and the velocity distribution, both for Maxwell molecules and hard spheres. Before that, the hydrodynamic profiles 



obtained from the simulations by using the two types of boundary conditions considered in Sec. Ill are presented 
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A. Hydrodynamic profiles 



As said in Sees. | and [fj], the Boltzmann equation for Maxwell molecules admits an exact solution characterized 
by Eqs. ([ri|)-([l9|). This solution applies to the bulk region, i.e., the region where boundary effects are negligible. 
Obviously, in a simulation with a finite size of the system, it is not possible to avoid boundary effects completely. 
On the other hand, one can expect that the "nonequilibrium" boundary conditions based on the BGK distribution, 
Eq. (^TJ), inhibit the boundary effects, as compared with the conventional "equilibrium" boundary conditions (|4C|). 
We have confirmed that this is indeed the case. As an illustrative example, let us consider A = 5 and a' = 0.92 for 
Maxwell molecules. This corresponds to 7bgk = 0.21, L = 4.68, Ul — 3.90, and 6l = —2.37. Figure || shows the 
temperature and velocity profiles as obtained by using the MB and BGK boundary conditions. It is apparent that the 
velocity slips and the temperature jumps at the walls are much larger in the former case than in the latter. Note that 
the maximum temperature is not exactly located in the layer adjacent to the lower plate but it is slightly shifted. It is 
interesting to remark that when plotting the temperature T as a function of the flow velocity u x , a parabolic curve is 
observed with both boundary conditions, as expected from Eqs. ([l8]) and (|l9|). A more evident proof of the advantage 
of the BGK conditions is shown in Fig. [s]. Since the pressure is a constant in the exact solution valid in the bulk, 
any deviation from p = const can be attributed to boundary effects. The pressure obtained with the BGK boundary 
conditions is practically constant, except near the upper plate, while the one obtained with the MB conditions is only 
nearly constant in a small region around y/L ~ 0.75. Finally, Fig. ^| shows the ratio a/a' between the actual (local) 
shear rate a measured in the simulations, cf. Eq. (|5^), and the imposed shear rate a' . The ratio is closer to 1 in the 
case of the BGK boundary conditions than in that of the MB conditions. In addition, in the former case the region 
where a is practically constant extends to layers closer to the lower plate. 

In summary, the above example illustrates that the BGK boundary conditions are much more efficient than the 
MB ones to measure transport properties in the bulk. Therefore, in what follows we will only consider the BGK 
conditions. In each case, we identify a bulk domain comprised between the layers y — yo and y — y\ where a ~ const, 
p ~ const, and 7 ~ const, and take averages of a, F, v F K , ^1,2, and $ over those layers. Typical values are yo/L ~ 0.2 
and yx/L ~ 0.8. 

B. Nonlinear transport coefficients 

In this subsection we compare the simulation results for Maxwell molecules and hard spheres obtained from the 
Monte Carlo simulations with the (universal) predictions of the Grad method and the BGK and ES kinetic models. 
As said in the Introduction, we are interested in situations where the Knudsen number is not small, so that nonlinear 
effects arc important. cj An interesting quantity in the nonlinear Couette flow is the parameter 7, which is related 
to the curvature of the temperature profile. Its shear-rate dependence is shown in Fig. |^. A remarkable feature is 
that the simulation data for both interactions seem to lie on a common curve. This indicates that, as predicted by 
the models, the transport properties are hardly sensitive to the interaction potential, provided that the quantities are 
conveniently scaled. While the kinetic models exhibit a good quantitative (in the case of the ES model) or qualitative 
(BGK model) agreement, the Grad method fails, except for small shear rates. Since 7(a) = a 2 F T] (a)/5F K (a), one can 
interpret 5j(a)/a 2 as an effective, shear-rate dependent Prandtl number (relative to the usual Pr). This quantity is 
bounded between 1 for small shear rates and 5/3 (in the BGK model) or 5/2 (in the ES model) in the limit of large 
shear rates. This is consistent with the fact that the BGK model underestimates the value of 7. 

The most important transport coefficient is the nonlinear viscosity represented by the function F r/ (a). According to 
Fig. |[ the three theories retain the qualitative trends of the simulation data, namely the decrease of F v with increasing 
a (shear thinning effect). In general, however, the kinetic models (especially the BGK model) have a better agreement 
than the Grad method. It is also interesting to remark that a slight influence of the interaction potential seems to 
exist, the shear thinning effect being a little bit more significant for hard spheres than for Maxwell molecules. The 
nonlinear thermal conductivity F K (a) is plotted in Fig. |9j. It is quite apparent tiat Grad's solution does not capture 
even the qualitative shape of F K , as was already noted in the case of hard disks.0 Again, the kinetic models present a 
good agreement, especially in the case of the ES model. 

Normal stress differences are characterized by the viscometric functions 5 , i,2(a)- These quantities are well described 
by the kinetic models, as shown in Figs. |o| and [ll]. At a quantitative level, the agreement is better in the case of 
Maxwell molecules. In fact, the viscometric functions, especially the second one, exhibit a certain influence of the 
potential. Although the functions ^1,2(0) were not evaluated from the Grad method in Refs. [To| and |l(j|, the friction 
function F^(a) was calculated. This quantity is plotted in Fig. [12], showing an agreement between the theories and 
the simulation data similar to the one found in Fig. || for the viscosity function. The last transport coefficient is 
the cross-coefficient $ defined by Eq. (^). The comparison with simulation results for this quantity is a stringent 
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test of the theories, since it is a generalization of a Burnett coefficient that measures the component of the heat flux 
orthogonal to the thermal gradient. Figure shows that, as happened with the thermal conductivity function F K (a), 
the Grad method gives a wrong prediction for the shear-rate dependence of <&(&). On the other hand, the kinetic 
models describe fairly well the nonlinear behavior of this function. In the case of the ES model, the agreement with 
the simulation data is practically perfect. 

C. Velocity distribution function 

Apart from the transport coefficients, the kinetic models provide the velocity distribution function. In the case 
of the BGK model, the solution is given by Eq. (|33|), while the ES distribution function can be found in Refs. || 
and ||. In order to assess their reliability, we have computed in the simulations the marginal distribution (p8|) at the 
layer y/L = 0.5. The ratio R(£ y ) = f{i y )/[-K^ 1/2 exp(-^)], where ip(£ y ) is defined by Eq. ©, is a measure of the 
departure of the distribution function from the local equilibrium. This quantity is plotted in Fig. [l4| for Maxwell 
molecules in the case of a reduced shear rate a — 0.636 and a reduced (local) thermal gradient e = —0.272. Both 
theories capture the main features of the actual distribution. While the BGK distribution exhibits a better agreement 
near the maximum (around £ y ~ —0.5), the ES distribution seems to describe better the region £ y > 1. The case of 
hard spheres is illustrated in Fig. |lj, which corresponds to a = 0.419 and e = —0.195. In this case, the ES model 
shows a superiority over the BGK model both near the maximum and for large positive velocities. 

V. CONCLUDING REMARKS 

This paper has dealt with the steady planar Couette flow in a dilute gas beyond the scope of the Navier-Stokes 
description. This nonlinear problem had been previously studied by means of kinetic theory tools, such as the Grad 
method, and the BGK and ES kinetic models. These theories predict momentum and heat fluxes characterized by five 
shear-rate dependent generalized transport coefficients: a viscosity function F v (a), Eq. (pp|), a thermal conductivity 
function F K (a), Eq. (pl|), two viscometric functions ^i^(a), Eqs. ( f22|) and (^3|), and a cross coefficient ^(a), Eq. (|25|). 
The main motivation of our study has been to perform DSMC simulations for Maxwell molecules and hard spheres 
in order to assess the reliability of the above theories in the non-Newtonian regime. Since we have been interested in 
the bulk properties, we have used "nonequilibrium" boundary conditions to inhibit the influence of finite-size effects. 

An important outcome of the simulation results is that, as predicted by the kinetic theories considered here, the 
shear-rate dependence of the transport coefficients is practically insensitive to the interaction potential when the 
quantities are properly nondimensionalized. In particular, the actual shear rate has been reduced with respect to an 
effective collision frequency defined from the Navier-Stokes shear viscosity coefficient. The simulation results show, 
however, that the second viscometric function presents a non negligible influence of the interaction model, so that the 
normal stress difference P zz — P yy is smaller for Maxwell molecules than for hard spheres. The comparison with the 
theoretical predictions shows that the kinetic models give a fairly good description of the five transport coefficients. 
On the other hand, the Grad method yields a shear viscosity in qualitative agreement with the simulations but it 
dramatically fails for the coefficients measuring the heat flux. This is basically due to the truncation scheme of the 
Grad method at the level of the heat flux. The physical idea behind a kinetic model is quite different, since it consists 
of replacing the true Boltzmann collision operator by a simple relaxation term but otherwise all the velocity moments 
are taken into account. As a consequence, while in the Grad method one has to solve a closed set of coupled differential 
equations for the moments, in the case of the kinetic model one gets the velocity distribution function and determines 
the fluxes from it. 

In the ES kinetic model the reference distribution function appearing in the collision operator is an anisotropic 
Gaussian parameterized by the pressure tensor. This allows the model to give the correct Prandtl number Pr = | 
but at the expense of complicating its mathematical structure. In the case of the BGK model, however, the reference 
distribution is that of local equilibrium but the model leads to Pr = 1. The agreement with simulation of the ES model 
is generally better than that of the BGK model, especially in the case of F K and $. In spite of this, it is fair to say that 
the performance of the BGK model is quite good, given its simplicity relative to that of the ES model. Finally, the 
results reported in this paper clearly shows the usefulness of kinetic models to analyze nonlinear transport phenomena 
in the Couette flow problem. This complements previous conclusions drawn from other nonlinear problems, such as 
the uniform shear flow and the Fourier flow. 
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APPENDIX: THEORETICAL EXPRESSIONS FOR THE TRANSPORT COEFFICIENTS 

In this Appendix we list the explicit shear-rate dependence of the dimensionless transport coefficients defined in 
Sec. |0], according to the Grad method, the BGK model, and the ES model. 

1. The Grad method 

From the Appendix of Ref. [l(], corrected in Ref. [l6|, one has 

w = i + ii + A (a r (A1) 

^ (a) = l-fJ + 3A( fl ) ' ^ 
1 - -26-a 2 

* ( "> = - 7 i + y + (1 - 8J)A(o) - («) 

where A(a) = \J^ + — |ff a4 ' viscometric functions are not evaluated in Refs. |l^ and |l6|, although the 

friction function is provided. It is given by 

F » = TTT2 v ( A4 ) 

1 + ^far + A(a) 



Note that F v (a) and F )Ji (a) become meaningless for a 2 > 25(vl057 4- 29)/432 ~ 3.56, while F K (a) and $(a) are 
unphysical for a 2 > 50/63 ~ 0.79. For small shear rates, the above transport coefficients behave as F v ~ 1 — ^a 2 , 



F K ~ 1 + §a 2 , $ ~ -|(1 - ±§a 2 ), and F„ ~ 1 - fa 2 . 

2. The BGK model 

The derivation of the transport coefficients from the BGK modefl implies the resummation of asymptotic series by 
means of the Borel method. As a consequence, the results are expressed in terms of the functions F r {x) defined by 
the recurrence relation F r (x) = [(d/ dx)x] r Fq(x) , where 

F (x) = - / dtte-^^K^t^/x 1 ^), (A5) 
x Jo 

Kq being the zeroth-order modified Bessel function. The curvature parameter 7(a) is given by the solution to the 
following implicit equation 

2 3Fi + 2F 2 

a = 7 p ' ( A6 ) 

where the functions F r are evaluated at x = 7. Next, the transport coefficients are expressed in terms of a and F r ( , y) 
as 

F v (a) = F , (A7) 
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. , F a 3Fi + 2F 2 . . 

F K (a) = — , (A8) 

* 1 (a) = -2F 1 3Fl+4F \ (A9) 
V ' 3Fi +2F 2 K ' 

*(o) = - [5F 2 + 2F 3 + a 2 (F 2 + 5F 3 + 8F 4 + 4F 5 )] . (All) 

For small shear rates, one gets F„ ~ 1 - fa 2 , F K ~ 1 - ^a 2 , ¥ x ~ -^(1 - ^ a 2 ), * 2 - f (1 - ^a 2 ), * - 
and F M ~ 1 - ^a 2 . 

3. The ES model 

In Refs. H and |j] the solution of the ES model is worked out keeping the Prandtl number Pr as a free parameter. 
Here we particularize the results to the correct value Pr = | . It is convenient to express the transport coefficients in 
terms of an auxiliary parameter (3, defined as the solution of the implicit equation 



fl2 = 4/3 [2f3{F 1 + 2F 2 ) - 3] 2 [3F + 2F 2 - 2£F\{F\ + 2F 2 )] 
9 F 2 [2/3(F 1 +F 2 )-3]+F 1 [2/3(F 1 +2F 2 )-3] 2 ' 



where now the functions F r are evaluated at x = p. The relationship between the curvature parameter 7 and (3 is 



7(o) = |/3[3-2/3(Fi + 2F 2 )]. (A13) 



The transport coefficients are 



In Eq. (|A18|) 



9Fn 

FJa) = — —2— - (A14) 

" V ; [2/3(Fi + 2F 2 )-3] 2 ' ^ ' 

F K {a) = T^F^a), (A15) 
57(a) 

m r \ - 12/3 3F 1 +4F 2 -2(3F 1 (F 1 +2F 2 ) 
l[a) ~ a 2 (3-2/?F 1 )[3-2/3F 1 (F 1 + 2F 2 ]' [ ' 

y 9 (a) = ^ El fA17l 

2W a 2 (3 -2/3F x )[3- 2(3F 1 (F 1 + 2F 2 y { 1 

$(a) = ic x {o 2 [Cf F 3 (Fi + 2F 2 ) - 9C 3 F 2 (F 2 + 2F 3 ) + 54C 2 F (F 2 + 4(F 3 + F 4 )) 

-108Ci(F 2 + 5F 3 + 4(2F 4 + F 5 ))] - ^(CiFoF - 6F 2 ) + 4C 2 F (F! + 2F 2 ) 

+4Ci[C 2 F Fi - 6(F 2 + 2F 3 )] - 24C 2 F 2 } . (A18) 



c ' s + (A19) 
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C 3 = Za 2 C(F£ + 4C 2 {2f3[C 1 (F 1 + 4F 2 ) + 3Fi] - 3d}. (A21) 
For small shear rates, one has F v ~ 1 - fa 2 , F K ~ 1 - ^a 2 , *! ~ (1 - ff a 2 ), * 2 ^ f (1 - ^a 2 ), * -£, 



and F M ~ 1 - 3a 2 . 
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FIG. 1. Plot of the marginal distribution function of particles reemitted by the walls. The solid line corresponds to the 
distribution measured in the simulation by applying the BGK boundary conditions (with a! = 1, eo = 0, and tL = —3.15), 
the dashed line is the theoretical BGK distribution (which is hardly distinguishable from the solid line), and the dotted line 
corresponds to the MB boundary conditions. 
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FIG. 2. Temperature and velocity profiles for hard spheres in the case L = 1.81, Ul = 3.17, and To — 6. The solid lines are 
results obtained from a DSMC simulation of the BGK equation, while the dashed lines are the theoretical BGK results. 
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FIG. 3. Marginal velocity distribution function for hard spheres at y = 0.5L in the case L = 1.81, Ul = 3.17, and To = 6. 
The solid line is obtained from a DSMC simulation of the BGK equation, while the dashed line (which is hardly distinguishable 
from the solid line) is the theoretical BGK distribution. 
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FIG. 4. Temperature and velocity profiles for Maxwell molecules, as obtained from DSMC simulations, in the case L — 4.68, 
Ul = 3.90, and To = 6. The solid lines refer to the results obtained from the BGK boundary conditions and the dashed lines 
refer to those obtained from the MB boundary conditions. 
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FIG. 5. Same as in Fig. 0, but for the pressure (which is measured in units of nfcsTi,). 
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FIG. 6. Same as in Fig. ^, but for the ratio a/a' between the actual shear rate, a, and the one imposed by the BGK 
boundary conditions, a' . 
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FIG. 7. Plot of the thermal curvature parameter 7 as a function of the shear rate. The symbols are simulation data for 
Maxwell molecules (circles) and for hard spheres (triangles), while the lines are the theoretical predictions given by the ES 
kinetic model (solid line), the BGK kinetic model (dashed line), and the Grad method (dotted line). 
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FIG. 8. Same as in Fig. (?], but for the viscosity function F v . 
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FIG. 9. Same as in Fig. M, but for the thermal conductivity function F K . 
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FIG. 10. Same as in Fig. M, but for the first viscometric function ^i, relative to the Burnett value 5'i(0). Note that the 
Grad prediction is not plotted. 
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FIG. 11. Same as in Fig. M, but for the second viscometric function $2, relative to the Burnett value ^2(0). Note that the 
Grad prediction is not plotted. 
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FIG. 12. Same as in Fig. H, but for the friction function _F M 
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FIG. 13. Same as in Fig. M, but for the cross coefficient $, relative to the Burnett value $(0). 



24 



R 



1.50 



1.25 



1.00 



0.75 



0.50 



0.25 -■ 



0.00 




-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 



FIG. 14. Marginal velocity distribution function, relative to the local equilibrium distribution, for Maxwell molecules at 
y = 0.5L in the case a = 0.636 and e — —0.272. The solid line is obtained from a DSMC simulation, while the dashed line is 
the theoretical 
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FIG. 15. Same as in Fig. [H| but for hard spheres in the case a = 0.419 and e = —0.195. 
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